%% process_fespace
function fespace = process_fespace(mesh, basis_type)
[fespace.node, fespace.elem] = process_refine(mesh, basis_type);
fespace.basis_type = basis_type;
fespace.Nb = size(fespace.node,1);
fespace.Nlb = size(fespace.elem,2);
end

%% process_refine
function [new_node, new_elem] = process_refine(mesh, basis_type)
switch basis_type
    case "P1"
        new_node = mesh.node;
        new_elem = mesh.elem;
    case "P2"
        midpoint = 0.5*(mesh.node(mesh.edge(:,1),:) + mesh.node(mesh.edge(:,2),:));
        new_node = [mesh.node; midpoint];
        e1 = abs(mesh.elem_ed)+mesh.Np;
        new_elem = [mesh.elem, e1];
    otherwise
        error("Invalid basis type.");
end
end